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Scope of Contract NAS2-415I 


Work under Contract NAS2-4151 started on February 1, 1967 
with the purpose T, to define and study one or more probabilistic 
models for the computation of . lifting rotor random dynamic 
loads and vibrations, which are suitable for the interpretation 
of wind tunnel and flight measurements of rotor loads and 
vibrations using power spectral density measuring techniques . u 
This report summarizes the results obtained through August, 1967 
when eight man-months of work had been expended. Work under 
subject contract, which extends through July, 1968 to a. total 
of twenty-one man-months, is being continued. 
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Concepts for a Theoretical and Experimental 
Study of Lifting Rotor Random Loads 
"and Vibrations 


by Kurt H. Hohenemser 
and Gopai H. Gaonkar 
Washington University, St . Louis, Missouri 


Abstract : 

After briefly discussing a number of lifting rotor conditions 
with random inputs, the present state of random process 
theory, applicable to lifting rotor problems is sketched. 
Possible theories of random blade flapping and random blade 
flap -bending are outlined and their limitations discussed. 

A plan for preliminary experiments to study random flapping 
motions of a see-saw rotor is developed. 
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1 . Introduction. The Significance of Lifting Rotor Random Loads 

As of now the problem of lifting rotor dynamic loads 
and vibrations has been treated almost exclusively on the 
basis of deterministic modeling. In the most advanced treatment 
the flow field In the vicinity of the rotor is computed from 
the systeta of vortices generated by the blades, and the dynamic 
response of the blades is obtained by computing the interactions 
between the flow field and the blades. In comparing computed 
and measured dynamic loads and vibrations, it. is necessary to 
assume stabilized flight conditions. Maneuver and gust dynamic 
loads and vibrations are treated a3 transients between stabilized 
conditions. 

This deterministic model of the generation of dynamic 
lifting rotor loads and vibrations is in many respects unsatis- 
factory, Even under the steady flow conditions for lifting 
rotor models in wind tunnels the blades perform sizable random 
motions indicating the gross turbulent character of the rotor 
inflow. In flight this local turbulence of the inflow is super- 
imposed to the atmospheric turbulence and to pilot control 
inputs which are both of a random nature. It would, therefore, 
seem appropriate to try to treat the problem of lifting rotor 
dynamic loads and vibrations with probabilistic modeling techniques. 
Such techniques have been successfully applied to the problem of 
airplane gust responses, see for example Ref. (l) and (2). 

Airplanes with known dynamic response characteristics are used 
to measure parts of the atmospheric turbulence spectra, and 
these spectra are then used to predict the gust responses of 
new prototype airplanes in the design stage. It is conceivable 
that an approach along these lines could be very useful for 
the prediction of lifting rotor dynamic loads and vibrations, 
though the problem Is here much more involved than for airplanes, 
because of the numerous significant degrees of freedom and the 
complexity of the random phenomena associated with the rotary 
wing aircraft . It is rather surprising that apparently no 
previous efforts in. this direction have been made although for 
some widely differing lifting rotor problems a description of 
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the underlying phenomena by random process is called for. 

The only indication that work in this area is planned was found 
in Ref. (3) which describes a test set-up suitable for measuring 
the response of a beam to lateral random loads. 

The stochastic methods developed to analyze airplane 
responses to atmospheric turbulence or to define atmospheric 
turbulence by airplane response measurements are not applicable 
to lifting rotors for two reasons: 

First, the coefficients of the equations of motion 
of lifting rotors are time variable rather than constant as 
for airplanes. Second, the Widely made assumption that the 
vertical gust velocities over the wing span can be approximated 
by a single stationary random process is not valid for rotary 
wings. The rotor blade will have to be subdivided in a number 
of elements whereby the vertical velocity components at' different 
elements will be represented by different non- stat lonary 
correlated random processes . The multiplicity and non- 
stationarity of the random processes associated with rotary 
wings invalidates the widely used relation according to which 
the ratio of output over input power spectral density equals 
the square of the system transfer function. Unlike airplanes, 
it -will be necessary to consider variations of the turbulence 
spectrum along the span of the blade. The problem then would be 
to derive Buch a distribution of the turbulence spectrum over 
the blade span from the measured blade motions, or in a' design 
situation, to derive the blade random vibrations from a given 
distribution of the turbulence spectrum over the blade span. 

For this reason neither the theoretical nor the experimental 
methods available for determining airplane responses to 
atmospheric turbulence can be adopted for lifting rotors, and 
new theoretical and experimental techniques must be developed. 

Six types of lifting rotor problems where the applica- 
tion of random process concepts appears to be promising will 
be briefly discussed. 
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1.1 Lifting Rotor Response to Atmospheric Turbulence 

In order to study the lifting rotor response to atmos- 
pheric turbulence one could think of testing a rotor model 
in a wind tunnel equipped with a gust generating system. 

Such a system has been installed In the NASA Langley Transonic 
Wind Tunnel and is described in Ref. (4). The system simulates 
a low amplitude continuous sinusoidal vertical gust. Under the 
usual assumption made for frozen wing aircraft - the vertical 
gust velocity can be represented by a single stationary random 
process - it is sufficient to measure the model response to 
sinusoidal vertical gusts over an adequate frequency range, 
since one needs only the system transfer function. For lifting 
rotors, however, for reasons explained before, the knowledge 
of the transfer function is not sufficient. The gust generator, 
in order to be applicable to lifting rotor models, must 'be 
capable of producing the entire gust spectrum with the proper 
amplitude ratios. No such device has apparently been designed 
or built anywhere as yet. Correspondingly, a theoretical 
computation of system transfer functions, as it is performed 
by the various sophisticated computer programs developed for 
lifting. rotors, is also not adequate for the problem of lifting 
rotor response to atmospheric turbulence, and considerably 
more elaborate programs will have to be developed for this 
purpose. 

1.2 . Stopping and Folding of Lifting Rotors in Flight 

The helicopter industry is vigorously working on 
problems of stopping and folding of lifting rotors in flight 
in an effort to develop air vehicles which - combine the 
hovering efficiency of the helicopter with the cruising 
efficiency of high performance airplanes i Wind tunnel tests 
have been conducted to study the stopping and folding process 
of lifting rotor models in a steady flow. Because of the , 
crucial effects of atmospheric turbulence these tests cannot 
be considered realistic. Even if the; transfer functions were 
to be determined analytically or experimentally in a wind tunnel 
with low amplitude sinusoidal gust simulation, the results 
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could not be applied to the problem (except when the rotor is 
stopped) because of the frequency shifts produced by the time 
varying coefficients in the equations of motion. The solution 
of the problem requires new methods of analysis, and in an 
experimental approach the model must be subjected to often ^ 
repeated stopping and starting conditions with a simulation of 


the full turbulence spectrum. 

1. 3 Lifting Rotor Hover ing Performance_ 

In spite of several decades of lifting rotor research 
and development, hovering performance prediction for new types 
of rotors is still quite poor. The widely accepted and contin- 
uously refined methods of performance prediction very often 
miss the mark by relatively large margins. A 10* overestimate 
of hovering performance is usually equivalent to 30 to 
overestimate of payload and can have very serious consequences . 

It is also known that certain types of atmospheric turbulence 
can considerably degrade the helicopter OGE hovering performance. 

Great progress has recently been made in the computation 
of the vortex wake system generated by a hovering lifting rotor. 
These computations, in agreement with test observations, have 
shown that the trailing vortices of one blade can get into the 
path of the subsequent blade whereby considerable drag increases 
and lift losses occur. The presence of atmospheric turbulence 
will amplify this tendency. A satisfactory definition and 
prediction of lifting rotor hovering performance should, 
therefore, include the effects of atmospheric turbulence and 
the associated random control inputs. 

1.4, . Unsteady Blade, Aerodynamics 

Another area of lifting rotor research where lately 

considerable progress has been made is that of unsteady blade 
aerodynamics, particularly in the region of partial blade 
stall. The unsteady stall phenomena are highly complex and 
very much dependent on the time history of blade angle of attack 
changes. Atmospheric turbulence will probably have a large 
effect on unsteady blade stall. Since the blade stall phenomena 
provide the forward speed limitation of most rotary wing 
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aircraft, theory and testing in this regime also seems to 
require the application of the concepts of random processes. 

1.5 Vertical Descent in the Vortex State 

Another flight condition where the lifting rotor 
loads are basically of a random nature, is the descent in 
the vortex state. Both rigid and flapping rotors show large 
thrust variations beginning at quite small vertical descent 
rates and becoming maximum at the fully developed vortex 
state inflow pattern. Such thrust variations have been 
described for example in Ref. (5). A theory using random 
process concepts would appear to be suitable for this flight 
condition. 1 

1.6 Transition from Hovering to Forward Flight and Vice Versa 

It is well known that most helicopters have a slow 
speed range corresponding to advance ratios between .05 and .10 
where dynamic blade loads and vibrations are much higher than 
in cruising flight. Take-off and landing maneuvers fall into 
this slow speed range and contribute to the high level of 
dynamic loads and vibrations. Flow through the rotor disk 
is Strongly non uniform and wake recirculation takes place. 

No adequate deterministic theory of this complex flow state 
has been established as yet. The transition phenomenon also 
appears to be a promising candidate for the application of 
random process theory. 

The preceding six cases where lifting rotor responses 
to random inputs are of importance are by no means exhaustive. 

It is also not at all obvious how to approach these problems. 
Their listing here has merely the purpose to provide some 
motivation to the search for suitable methods of treating lifting 
rotor random loads and vibrations. 



2. Survey of Random Process Theory Related to Lifting Rotor 
Problems 

Random process theory so far has been predominantly 
applied in communication and automatic control engineering. 

Its application to mechanical systems and structures is still 
quite limited. A number of good textbooks on random process 
theory and measurements have recently been published, for example 
Ref. (6) and (7), and some modern textbooks on mechanical 
vibrations contain a chapter on the response of structures to 
random inputs, see for example Ref. (8). Most of the 
literature deals with constant parameter linear systems subjected 
to a stationary random input. Very little work exists on non 
stationary random inputs as they occur in lifting rotors, 
and no previous work was found on the computation of responses 
of time varying linear systems to 'random inputs which is 
required in the random flapping and flap-bending analysis of 
lifting rotor blades. The response of complex structures with 
many degrees of freedom subjected to multiple correlated random 
inputs has only very recently been treated in a few papers, 
see for example Ref.‘ ( 9 ) and (10), and the best approach to 
this problem which must be solved for the lifting rotor, is by 
no means obvious. This section gives a brief survey of random 
process theory as related to lifting rotor problems. 

2.1 General Remarks on Random Processes 

A random process is defined as an ensemble of random 
sample functions x(t) which can be described by a variety of 
probabilistic measures. The first' order probability p(x-| ) 
gives the distribution of Xj = x(t^) over the ensemble at the 
time t-j . For the sake of simplicity we will assume here that 
the expectation or mean value of x(t) is always zero: 

oo 

fl x = E[x(t)] * j xp(x)dx e 0 2.1.1 

— oo 

The second order probability p(x^,x 2 ) gives the joint distri- 
bution of x.j » x(t-j) and x 2 = x(t 2 ) at the times tj and t 2 * 

The expectation of the product x^x 2 is called the autocorrelation 
function R x (t 1 ,t 2 ) 
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oo 

R x (t 1 ,t 2 ) = E [ Xl x 2 ] = jfx 1 z 2 p(x 1 ,x 2 )dx 1 dx 2 

— 05 

and, is in general a function of both and t 2 . 

Higher order probabilities p(x^ ,x 2 ,x 3 , . . . ) and higher 
moments E [ x.|X 2 x 3 . . .] can be defined and are necessary for a 
complete description of a random process. Usually one limits 
the description to the first and second order probabilities. 

For the special case of a normal or Gaussian process, such a 
description is complete. 

In the case of two correlated random processes with 
sample functions x(t) and y(t) there exists a second order joint 
probability p(x<j>y 2 )', where x-j » x(t^); y 2 « yCtg)* 
expectation of the product x^y 2 is called the cross-correlation 

function 

t x 1 y 2 3 “ 

Vie now introduce the complex sample Fourier transforms of the 
sample functions x(t),y(t) 


O0 

Jj x 1 y 2 p(x 1 ,y 2 )dx 1 dy 2 
— 00 


2 . 1.3 


X(f) = J x(t)e' 12lrft dt 

-OO 

Y( f ) = J y(t)e" i2irf dt 


with their inverse 


x(t) = f X(f)e i2irft df 

— CO 
OO 

y(t) - J Y(f)e i2irft df 


2 . 1.4 


2 . 1.5 
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The probabilistic measure in the frequency domain 
which is equivalent to R (t-^tg) in the time domain is the. 
expectation of the product of complex sample Fourier transforms 
E [ X*(f -j)Y(f 2 ) ] which is called the double frequency cross- 
spectral density . It can be shown that it is related to the 
cross-correlation function by 

00 

s (f v f 2 ) - E [x*(f 1 )Y(f 2 ) ] = JJ R xy (t 1 ,fc 2 )e i 2 ,r ( f 1 t r f 2 t 2ht 1 dt 2 

J -oo 

2 . 1.6 


The inverse relation between cross-correlation function and 
double frequency cross-spectral density is given by 



jj S xy (f 1 ,f 2 )e- i2lr(f 1 t T f 2 t 2 ) df 1 df 2 

— OO 


2 . 1,7 


The corresponding relations between auto-correlation function 
and power spectral density are obtained by replacing by R x , 

S xy by S x and Y(f 2 ) by X(f 2 ). 

Random processes are called weakly stationary if the 
first order probabilities are independent of time and if auto- 
correlation and cross-correlation functions depend only on the 
time difference r = t 2 ~t.j and no ^ on In ^his case 

one can prove the inequality 

|R xy ( r) J 2 < ^(0)^(0) 2.1.8 

which allows to define the normalized cross-correlation function 


P 

xy 


( r ) « 



( r ) 


Vr x (°)V°> 


R xy< r > 

<*x. 


2 . 1.9 
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a are the standard deviations of the time 

y 

For 


where <r x and 
independent probability distributions p(x) and p(y). 


x = y one obtains 


lR x (r)| < |H X (0)1 = <jr x 


2.1.10 


r a 


If in the time domain R^, (tj ,t’ 2 ) depends only on 
and not on t^+tg, it can be shown that in the frequency 
domain S (f«,f P ) also reduces to a function of a single 
frequency given by 


Vh 


oo 

S xy (f) = E [x*(f)Y(f) ] = J R xy ( T )e“ i27rfr dr 

“00 


2 , 1.11 


For x » y this reduces to 

S x (f) = E [x*(f)X(f)] 


2.1.12 


In words: The power spectral density of a stationary random 
process is the expectation of the product of a sample Fourier 
transform with its conjugate complex value. 

If the expected values taken over the ensemble of 
random functions are equal to the corresponding time averages 
taken over a sample function, the stationary process is called 
ergodic. In this case 


R xy* = E t x 1 y 2^ 


lim 
T— oo 2T 


j . 

i j x (t)y(t+ r )dt 
-T 


2.1,13 


While the double frequency' power spectral densities of non- 
stationary random processes cannot be measured unless an entire 
ensemble of random functions is available, S xy ( f ) and S x (f) 
for stationary ergodic processes are measurable, if only a 
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single sample function is given. In. order to obtain S (f) 

i 

one merely has to send the signal from a random sample 
function through a narrow band pass filter with the center 
frequency f and band width Af and square and average it over an 
adequately long time period and divide the result by Af . 

For a ergodic process - again assuming zero mean 
value of the random function - the standard deviation a 
of the probability distribution can be expressed either by the 
value of the autocorrelation function for r » 0 or by the 
integral over the power spectral density: 



lim Ji__ 
T-oo 2T 



oo 

= R x (o) = J s x (f)df 

—do 

2.1.14 


Thus, for an ergodic process, the standard deviation of the 
probability distribution over the ensemble of random functions 
at any given time t equals the root mean square value taken over 
any sample function and can be obtained also by integrating 
over the power spectral density distribution. For an ergodic 
random process the time rate of the random function is uncorre- 
lated to the function value at the same time. The standard 
deviation of the time rate is given by either one of the 
4 expressions: 

T 

OO - oo 

m J x 2 p(x)dx = J x 2 (l>)clt = -R x (0) = J (2TTf) 2 S x (f)df 

— OO _rj» —00 

'2.1.15 


An important concept applicable to the fatigue of structures 
under random loads is. the expected number of crossings of a 


certain level per unit time. For a normal random process 
one can show that the expected value of positive crossings 
per unit time of the level a is given by 





2.1.16 
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Making use of the two preceding equations the expected value 
of positive zero crossings (a = 0) per unit time is 



_2 

2t r 


J(2irf) 2 S x (f)df 



2.1.17 


For a narrow band process almost every positive zero crossing 
leads to a full cycle , so that this equation, when applied to . 
narrow band processes, also gives the expected value of cycles 
per unit time which is important in fatigue considerations. 

2.2 Response of a Constant Parameter Linear System to a 
Stationary Random Input 

A linear dynamic system with constant parameters can be 
described by the unit impulse response function h( t ) . For a 
general input x(t) the system response is then given by the 
convolution integral 

CO 

y(t) as f h( T)x(t-T )dT 2.2.1 

<T 

Taking the Fourier transform on both sides of this equation 
results in the relation 

Y(f) * H(f)X(f) 2.2.2 


where 

H(f) = j h{ r )e" i27rfT dr 2.2.3 

—oo 

is the complex frequency response function of the system. 

For stationary ergodic random inputs one determines the response 
auto-correlation function by taking the expected value of the 
product y(t)y(t+r). With the help of the convolution integral 
this results in the following relation between input and response 

■ i ■ . " ‘ . ' 

autocorrelation functions 
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R y ( T ) = jj"h( « )h( V )R X { T+f -»7 )d$ d»j 2.2.4 

0 

The corresponding relation in the frequency domain is 

S y (f) = I H(f) | 2 S x (f) 2.2.5 

All terms in this simple equation are real. The response 
standard deviation is given by 

CD 

<T y 2 = J I H(f) | 2 S x (f)df 2.2.6 

— 00 

and o^. is again equal to the root mean square value taken 
over a sufficiently long time period of a sample response 
function. 

The relation according to which the response power 
spectral density is equal to the input power spectral density 
multiplied by the square of the complex frequency response 
function is very widely used in many applications of random 
process theory. However, its validity is strictly limited 
to the case of a constant parameter linear system with a 
single stationary ergodic random input. 

2 . 3 Response of a Constant Parameter Linear System to a 
Non Stationary Random Input 

A general solution of this problem has been developed 
by Caughey in Ref . (ll) in form of complex double and triple 
integrals. The practical problem is that, except in special 
cases, the complex double frequency power spectra associated 
with non stationary random processes, cannot be measured with 
a single sample function, since it takes the entire ensemble 
of functions to define the variation of its probabilistic 
measures with time. There are, however, two special cases 
where the power spectral densities associated with non stationary 
random processes can be defined and measured when only a single 
sample function is available. In both cases the non stationarity 
is due t© a deterministic time trend. Such phenomena have been 
treated in Ref. (12). 
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In the first case the time trend is slow as compared 
to the random instantaneous fluctuations, so that the pro- 
babilistic measures can be determined with a single sample 
function by selecting an averaging time which is sufficiently 
long to obtain statistically meaningful values for the auto- 
correlation function or power spectral density but which is 
sufficiently short to exclude substantial effects of the slow 
time trend. Random processes which allow such measuring proce- 
dures are called locally stationary. 

In the second case the deterministic time trends are 
periodic so that time averages of autocorrelation function and 
power spectral density can be defined. Since this case is 
important in lifting rotor theory it will be discussed here in 
more detail. We assume that x(t) is a sample function of a 
stationary ergodic random process and we would like to know the 
average autocorrelation function and the power spectral -density 
of the modulated random sample function y(t) defined by 

y(t) = x(t.)cos27rFt 2.3.1 

The autocorrelation function is: 

Ryft^tg) = E [x(t 1 )x(t 2 )] COS 217^00827^2 
= (1/2)R ( r )(cos27rPT+cos2irP(t 1 +t 2 )) 2.3.2 

where r = t^-t.] • 

Averaging over a sufficiently long time period, the last term 
will disappear and one obtains for the time averaged auto- 
correlation function 

R„( r ) ■ ( l/2)R ( T )cos2irFT 

y x 


2.3.3 
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The average autocorrelation function of the modulated stationary 
random function x(t) is obtained by multiplying R { r ) with 
(1/2)cos2ttP . 

The corresponding power spectral density can be either 
obtained by using the relation 

oo 

S y (f) = f R x ( t )e" l2lrfT dt 2.3.4 

— OO 

or by taking the Fourier transform Y(f) of y(t) and inserting 
it in the definition 

S y (f) = E [y*(f)y(f)] 2.3.5 

The result is in either case 

s y (f) = (1/4) [ S x (f-F)+S x (f+F)] 2.3.6 


This power spectral density would be measured if a sample 
function y(t) were analyzed in the usual way by sending the 
signal through a band pass filter with a band width Af which 
is small as compared to F and by squaring and time averaging 
over a time which is large as compared to l/F. If the response 
power spectral density of the constant parameter linear system 
with the forcing function y(t) were determined in the same way, 
it would be given by |H(f) | 2 S^(f), where H(f) is again the 
complex frequency response function of the system. 

2 .4 Response of Time Variable Linear Systems to Random Inputs 
A linear system with time varying parameters - like a 
lifting rotor flapping blade - can be described by the unit 
impulse response function, which now is not only dependent 
on the time r elapsed since applying the unit impulse, but 
also on the time t at which the impulse had been applied: h( r,t). 
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In general the computation of h( T,t) is not possible in closed 
form and numerical methods are required to find approximate 
solutions. For a general input x(t) the system response is 
again given by the convolution integral 


y(t) = f h( r r }d 2.4.1 

<T 

The frequency response function is also time variable and 
determined by 


H(f,t) = J h( T ,t) e i2irfr dT 2.4.2 

— oo 

Taking the expected value of the product y(t^)y(t 2 ) one obtains 
with the help of the convolution integral the following relation 
between input and response autocorrelation function: 


R y (t r ,t 2 ) = 


Jj* h{ 4j»t-j)h{ fj , 1 2 ) Rjj. ( t -j - ^ , t 2 - >7 )d£ d tj 2.4.3 


The corresponding relation in the frequency domain is 


00 

s y (r v r 2 ) - JJ s x ( f f v )j*( $ n , v-t 2 )^ 2.4.4 

—oo 


where 

00 

J(f,f) = J H( { ,t)e i2lrft dt 2 . 4.5 

—do 

Apparently no- attempt has been made as yet to solve these 
equations for an actual problem. A solution would be extremely 
difficult, since already the unit impulsq response h( f ,t) 
is not given analytically but must be approximated by numerical 



methods. The solution then requires triple infinite integrals 
over non analytical complex functions. Furthermore, as men- 
tioned before, the input and output double frequency power spectra 
cannot be measured with a single sample function. However, 
if the variable parameter system is electronically simulated 
the input-output relation for power spectral densities can be 
experimentally observed under certain conditions. 

2.5 Response of Linear Systems to Multiple Random Inputs 

The application of random process theory to complex 
structures with multiple random inputs - even assuming 
stationarity - has not as yet progressed to a point where a 
good understanding of the various possible and necessary 
approximations has been achieved. The problem of random 
blade flap bending is therefore presented in Section 4 in 
considerable detail. A substantial simplification is possible 
if, as usually assumed, the damping of the structure is small 
so that normal mode cross damping terms can be neglected. 
Unfortunately, lifting rotor blades have large aerodynamic 
damping at least in the low frequency inodes and not all cross 
damping terms are negligible. In this case, it may be simpler 
to transform the second order differential equations of the 
structure into an equivalent system of first order equations 
with respect to time* The normal modes are then complex and 
more difficult to handle, but the frequency response function 
is much simpler than for a second order system* Section 4 
includes the general solutions both for a second order repre- 
sentation with real normal modes and for a first order repre- 
sentation with complex normal modes. 
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3 . Random Blade Flapping 

For sufficiently slow excitation the lifting rotor blade 
will respond mainly in its fundamental mode, that is in the 
blade flapping mode'.; For biadet attached to the hub with 
flapping hinges, this is a rigid blade mode. For blades without 
flapping hinges the mode involves flap bending, however, in 
first approximation one can substitute rigid flapping about an 
off-set hinge. In the following a very simple analytical 
blade model is assumed which is adequate to discuss the basic 
problems of a random loads and vibration analysis. 

3 . 1 Assumptions 

It is assumed that the atmospheric turbulence affecting 
lifting rotor blade flapping consists of vertical velocity 
fluctuations which are uniform over the rotor span and which 
have a wave length which is large as compared to the effective 
rotor diameter. This type is called one -dimensional isotropic 
turbulence and is frequently used for airplane turbulence 
response calculations, see for example p. 21 of Ref. (l) . 

In a first order blade flapping theory valid for 
small flapping angles 0 , zero flapping hinge off-set and for 
moderate advance ratios m > one obtains the following equation 
of motion ( see for example Ref. (13), eqn. (15)) 

j8 4* ^-(1+ fx sint )0 + ( 1+ — ^ — cost)0 = — g 7 — (1+ ^ sint) 3.1*1 


The equation has been linearized both with respect to the 
flapping angle 0 and with respect to the rotor advance ratio m . 
The time unit has been selected so that the rotor angular velocity 
n is one. The only blade parameter in the equation is the 
non-dimensional blade inertia number y which for practical 
rotors has values between 2 and 10. Actually, the blade angle 
of attack at should vary both with radius and with azimuth 
angle. However, in this simple mathematical model an average 
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value a has been assumed which is independent of radius and 
blade azimuth angle and which is directly proportional to the 
vertical gust velocity. Using ,7R as effective rotor radius, 
the relation- between vertical gust velocity w and a is: 


Thus, the average angle of attack a directly reflects the 
atmospheric turbulence and is under the usual assumptions a 
stationary random function. 

The stationary random process a , according to the 

blade flapping differential equation, is modulated by the 

8 ■ 

deterministic time function (1+ -j v sint) and this non- stationary 
random input acts on a linear system with periodically varying 
parameters. In a more refined analysis considering the varia- 
bility of the angle of attack with radius and with blade 
azimuth, a would be replaced by the product of w with a deter- 
ministic time function. 

From the physical aspects of the problem it is clear, 
that the assumption of a locally stationary process, discussed 
in Section 2.3 is not satisfied for 0 or asint. The time trends 
from sint and cost are not slow as compared to the gust velocity 
fluctuations. 

3.2 Small Advance Ratio 

For small advance ratio m the periodicity of the para- 
meters in the blade flapping equation can be neglected. 

In this case the average power spectral density of the random 
input can be easily computed, since the modulated angle of 
attack a sint has the average power spectral density 

1/4 [s ff (f- (f+ ^)] 3.2.1 


see Section 2.3 (the relation given in Section 2.3 for the 
modulating function cos27rFt is also valid for sin27rFt). 
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If H(f) is the complex frequency response function for the blade, 
which, for constant parameters, has the value 

H(f) 5 —! 

-(2irf) +12irf^ +1 3.2.2 

one obtains for the average power spectral density of the 
flapping angle 0 the value: 

% (f> 


In the absence of a solution for periodically varying para- 
meters it is not possible to quantitatively establish the 
advance ratio range m for which the above solution is a good 
approximation. It would appear, however, that an advance 
ratio of v ~ .1 should show a negligible effect of the 
M -terms in the factors of 0 and 0 . It is noteworthy that 
this advance ratio range includes the region of . transition 
from hovering to forward flight and vice versa, during which 
considerable random loads and vibrations are encountered. 

3.3 Moderate Advance Ratio 

Inasmuch as up to now the only non stationary random 
processes which are analytically manageable are stationary 
processes which are modulated by a time function (such pro- 
cesses are treated in Ref. (12)), one might try to find a 
solution to the blade flapping equation 3.1.1 by assuming 
that the flapping angle 0 represents a stationary random 
process multiplied by a periodic time function. Since the 
input random process is of this type it may not be unreasonable 
to expect that at least to a certain approximation the response 
random process is of the same type . When time averaging the 
autocorrelation function of such a process one finds that the 
result depends only on T «= not on or same as for 


H(f) 



2 ir 


b+ s <? (f+ 


M 


3.2.3 
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stationary random processes. In the frequency domain, time 
averaging amounts to multiplying the double frequency power 
spectral density with the delta function a{f 2 ~f«|)> thus 
reducing it to a single frequency power spectral density. 

This is shown for the general case of a periodic modulating 
function in Appendix A. Here we will consider the special 
case 

y( t) as x(t)cos27rFt 3,3.1 

treated in Section 2.3. x(t) is again a sample function of a 
stationary random process. According to eqn. 2,1.6 the power 
spectral density for y(t) is: 

S y (f 1 ,f 2 ) = E[y*(f|)Y(f 2 )] 3.3.2 


Applying the postulated rule, that time averaging corresponds 
in the frequency domain to multiplication of the double 
frequency power spectral density with one obtains 

for the power spectral density of the random process with time 
averaged autocorrelation function 




3.3.3 

or 

S y (f) = E [Y*(f)Y(f)] 

3.3.4 

Since 

Y(f) = l/2(X(f+F)+X(f-F)) : and 

3.3.5 

since 

E j" X*( f+F)X( f-F)J =0 

3.3.6 

one obtains 



S y (f) = 1/4 [s x (f+F)+S x (f-F)] 


3.3.7 
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in agreement with eqn . 2.3.6* In particular, we have for the 

i 

power spectral density of the time averaged random process the 
rule, that the sample Fourier transforms are uncorrelated 
for different frequencies: 

E JVu^UpJ =0 for f i * fg 3.3.8 

same as for stationary random processes. 

Let us now assume that the random process with sample 
function 0(t) from eqn. 3.1.1 is of the type, where the 
averaging rule applied. We first write eqn. 3.1.1 in the form 

(Sf+fc-i+a^ slnt) ^ +(c 2 +a 2 cost) =» ( c 3 +a 3 clnt)a 3.3.9 

This equation is Fourier transformed by considering the following 
function pairs 

0(t) 

cost /j(t) 

0(t) 

sint 0(t) 

?(t) 

One then obtains 

{ *t 

(fy 2?) 

-(V k )B(f 1 + ^) } +c 2 B(f 1 ) + ^ {B(f r Ii)+B(f 1+ y\ 

o 3 A(f 1 )+ 2?) +S ( f 1+ 1?)} 


— * B(f) 

— * (1/2) | B(f- ^i-)+B(f+ 2^“)^ 

<— + 27rifB( f ) 3.3.10 

— )b( f - 2j->B( f +dr>} 

«— . -(2irf) 2 B(f) 


3.3.11 



- 22 - 


The conjugate complex equation with f 2 instead of f-j reads: 
-(2irr 2 ) Z B*lr 2 )-C)2wi.r z $(£ 2 )+^Tr {(fg- ~) 


-(f 2 + ^)B*(f 2 + | ¥ )}+c 2 B*(f 2 ) + -| { B*( fg- ^)+B*{f 2 + ^r)} 

C3 5* (f2 ) + !§ {a*( V ^) + r ( f 2+ Ip)} 3.3.12 

Multiplying the last two equations with each other and taking 
the expectation of the product leads to an equation for the 
double frequency spectral density (f-pf^). Multiplying this 
equation further with 5(f 2 -f 1 ) reduces the double frequency 
spectra to single frequency spectra corresponding to the time 
averaged autocorrelation function. In this final equation 
the imaginary terms cancel each other: 


Sp (f) | (2Trf) 4 +(o 1 2TTf) 2 -2c 2 (2irf) 2 +o 2 S } 

+% (f+ {-4-(2ir(f+ --j- 1 2lr(f+ £? )+ ~T~ I 

a 2 

(f - I?) { -i- (27r(f_ h? ' + 


a T a 2 


lp)+ -f-} 


= C 3 2 S ff (f)+ { S« (f- 5f) + S a ( f+ ^} 


3.3,13 


We have here a functional equation between S^r (f) and (f) 

which can be approximated by a system of linear equations 

valid at the frequency points 0,. Af , 2Af . . . By inverting 

this system of linear equations (f) can be obtained from 

Sg- (f). The details of the computation are presented in Appendix B. 

Again, in the absence of a rigorous solution, of the 
problem it is not possible to establish the advance ratio 
range for which eqn. 3.3.13 is useful. An attempt has been 
made to compare the result of eqn. 3.3.13 for a number of 
cases with the results of simulator studies. As of now, 
the comparison is not conclusive. 
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Numerical computations of Sp (f) in Appendix B, 
assuming a blade inertia number of 7 « 4 and an advance ratio 
of v *= .3 have shown, that for typical power spectral density 
distributions S s (f) equation 3.3.13 gives a blade flapping 
power spectral density Sp (f) which is approximately identical 
to that obtained with equation 3.2.3. It would then appear, 
that the damping and stiffness variation in eqn. 3.1.1 is 
negligible up to an advance ratio of u = ,3. The reason for 
this result is that for 7 « 4 and m « .3 the numerical values 

p p 

for a-j , a a 2 and ~ are two orders of magnitude smaller than 
the: numerical values for and c^. 
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4. Random Flap Bending 

Though ultimately the problem of random blade flap 
bending will have to be solved for non stationary random 
inputs including the periodic parameter variations of the. 
system, the analysis in this section ignores the non stationary 
character of the input and the time variability of the para- 
meters. Nevertheless, the many degrees of freedom and the 
multiple correlated random inputs require a rather complex 
analysis. The usual variability of stiffness and mass over the 
blade length make a rigorous solution of the dynamic problem 
not fea sable' and numerical approximations must be developed. 

The method selected here makes use of influence coefficient 
and stiffness matrices derived from a certain discretization 
of the blade . 

Three types of analysis are outlined. For the first 
type the system response is computed directly with the general 
admittance matrix. For the second type the response is com- 
puted with the help of the real normal modes of the undamped 
system. Under certain conditions normal mode cross damping 
terms are negligible and the admittance matrix is diagonalized. 

If these conditions are not satisfied a third type of analysis 
is possible where the second order differential equations with 
respect to time are transformed into a first order system. 

The eigenfunctions of this system are complex, but the admittance 
matrix can be diagonalized for this system without any limita- 
tions. 

' The discrete blade element considered in the present 
analysis is a simple one - a rigid element in which the bending 
stiffness is simulated by introducing torsional springs at the 
ends. The matrix formulation of the problem is not much mbre 
complicated if, instead, elastic elements with constant or 
linearly varying strains are used as suggested for example in 
Ref. (14), 
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4.1 Derivation of Blade Stiffness Matrix 

In a rotating reference system the equation for the 
vertical deflections y (positive up) of a rotating blade In . 
hovering, using qua si static aerodynamics and neglecting coupling 
with chordwise bending and torsion reads (see for example 
Ref. (15)) 

(Ely M ) -{Ty» ) Vmy+acp |^y - aca pii^- 4.1.1 


Instead of solving this differential equation by 

numerical methods the blade is subdivided into a number of 

rigid elements connected by torsion springs in the way shown 

in Fig. 1. Mass, aerodynamic damping, and Imposed external 

forces are assumed to be concentrated at the element midpoint* 

In order to obtain the influence coefficient or stiffness 

matrix of the discretized blade, we assume forces acting 

\ n 

at the element midpoints and compute the vertical deflections 
y n at these points from the equilibrium equations. 

The horizontal force equilibrium equations 

T .+m x, =* 0 , n '== 1,2...m 4.1.2 

n+i n n n 


are uncoupled from the remaining blade equilibrium equations 
and can be summed to obtain the tension force T n : 



m 


E 

p = rtf-1 


2 

U m x 
P P 


4.1.3 


m designates the furthest out blade element at the blade tip. 

The moment equilibrim about the left endpoint of 
element n is expressed by the equation 


M h- M n + 1 + 


m ft 1 _ 

n+inn 


S nf1 1 n + 


F — 
n 2 


4.1.4 
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DIMENSIONS AND LOADING FOR ANY TWO ADJOINING 
ELEMENTS WITH LUMPED PARAMETERS OF MASS, OAMPING 
AND STIFFNESS. 



FORCE SYSTEMS ON A TYPICAL ELEMENT «. 
(DYNAMIC FORCES ARE SHOWN IN BRACKETS) 


FIGURE 1 DISCRETE ELEMENT IDEALIZATION FOR BLADE 
FLAP - BiNOINO 
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or, introducing the shearforce 


m 

S n+1 ~ £ 

. j = n+1 J 


4.1.5 


by the equation 


m 


M n- M n + 1 +T n*Vn = E F j +F n IT 

j = n+1 


4.1.6 


There are m equations of this type 


Introducing the spring constant of the torsion spring at the point 

1 . 


X n+~1I by 


El 


K 


n+1 


n+1 “ x n+r x n 


4.1.7 


the torsion spring equation reads 


M 


e e « 


n+ 1 n K 


n+1 

n+1 


4.1.8 


For the cantilever blade there are m equations of this type. 
For a hinged blade the first torsion spring equation is 
replaced by the boundary condition M-j = 0. 

Finally the relation between slopes e n and vertical deflection 

y is obtained; 

17 n 


y n+1 ~ y n 


1„ 
ft rc 

0 n 2 


+ e 


n+V 


'n+1 


4.1.9 
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There are m equations of this type, the first one having the 
form 


*= e. 


y i - ^ 

Thus the 3m unknowns y 


T ' ' w nv' ~ 1 

defined by the 3m equations 4.1.6, 4.1.8 and 4.1.9* 


. . ,M are uniquely 

I m 


0 1 . . . 0 , 

.ii t rn 

Since the 

last set of equations is independent of the first two sets. 


0i 


e. 


are calculated from equations 4.1,6 and 


1 # * * ~nr 

4,1.8, while eqn. 4.1.9 is required for the calculation of 

the influence coefficients. Once the influence coefficient 
* / 
matrix [w] with the elements 

(j • 

w 


has been computed, the stiffness natrix [k] is obtained by 
inversion, since 


M M = W 


4.1.10 


4.2 Response Analysis with General Admittance Matrix 

The equations of vertical motion for an assemblage 
of discrete blade elements have in matrix notation the form 

[m] jyj + [c] { y } + [k] jy} = j f| 4.2.1 

In subscript notation these equations read 

m k y k + °k y k + £ k k/j = f k> k “ 

[m] , [c] and [k] are the mass, damping and stiffness 

matrices respectively. For the blade flap bending problem 
mass and damping matrices are diagonal, however, the theory 
is developed here for the general case of non diagonal mass 
and damping matrices. Various aspects of the response analysis 
of linear systems under random excitations are presented in 
References (8), (9)> (10), (l6), (17)* The common features 
in all these investigations are the following: 
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i) discretize the continuum to obtain a tractable spectral 
representation of the exciting system of random loads 
in different elements by different correlated random 
processes * 

ii) under the hypothesis of weakly stationary and ergodic 
behaviour of the random process determine the output 
power spectral densities from admittance matrix and input 
power spectral densities. 

The complex admittance matrix for the system defined by 
equation 4.2 .1 is 


|h(o>)J = w^[m] + iw + [kj j~ ^ 4.2.2 

This expression replaces the frequency response function given 
in eqn. 3.2.2 for blade flapping, w has been substituted 
for 2?rf. The complex natural frequencies are obtained as the 
roots of the characteristic equation 


|A 2 [»] + X [c] + [k] | . 0 


4.2.3 


Under weakly stationary and ergodic assumptions, the 
input power spectral density and the admittance matrix 
completely describe the output power spectral density. The 
power spectral densities now occur in form of matrices. 

For example the matrix with elements Sj k (w) is 

defined by 


3 jk<"> - e[pj*(-)P K (-)] 


4.2.4 


where Fj(«) is the sample Fourier transform of the input 
random force acting on blade element j. The deflection or 
output power spectral density matrix £s^(«)| with elements 
is determined by 


[s^")] = [h*(«)J [S f (u,)j [H(-)j 


4.2.5 
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This relation replaces eqn. 2-2.5 for a single degree of 
freedom system for which the output power spectral density 
equals the input power spectral density multiplied by the 
square of the frequency response function. 

The autocorrelation functions of the single degree 


of freedom system is now also replaced by a matrix*^ For the 
input forces we have the cross-correlation matrix j^R (^J 
with the elements Rj k (r) defined by 


R 5k< T > - E [fj^yt+r)] 


4.2.6 


For T= 0 one obtains the force cross-correlation matrix 
between blade elements 

oo 

fR f (0)l = J [S f (<ojJa« 4.2 


and a similar expression for the deflection cross-correlation 
matrix 

^(0)] = J [S y (wJ]dw 4.2.8 

This expression replaces eqn. 2.2.6 for a single degree of 
freedom. The diagonal terms of |^R y (0)J represent the 
mean square deflections of the various blade elements and 
the off-diagonal terms represent the deflection cross 
correlations between different blade elements. 

In principle, given the input power spectral density 
matrix £s f (w)j and the complex admittance matrix [h(«)] , 

the output power spectral density matrix [s'M can be 
computed from eqn. 4.2.5.' It is noteworthy that this analysis 
can be performed without making any assumptions with respect 
to the damping matrix. But a major drawback is the manipulation 
of the admittance matrix containing complex., elements. These 
difficulties are avoided in a modal analysis for which the 
admittance matrix is dlagonalised . 
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4.3 Real Normal Mode Analysis 

For many structural problems it is adequate to perform a 
normal mode analysis with the modes of the undamped system. If 
the damping is distributed in a certain way, or if normal mode 
cross damping terms can be neglected, one obtains in such a normal 
mode analysis uncoupled differential equations for the normal 
coordinates. Eqn. 4.2.1 can be easily transformed in such a way 
that the factor [m] in the first term is replaced by the unit 
matrix ex] • With the definitions 


|f} = [m] {?} 


if i - or i it 

[5] . [.]■*[.] on h- [■]'*[*■] [»n. 

and omitting the bars over all symbols Eqn. 4.2,1 is obtained 
in the form 

{y} + [o]{y} + [k] {yj = {f} 4.3.1 

where w and (X) are symmetrical matrices. The normal modes 
of the undamped system are obtained from 

2 


[« ] iy ] - o ] pi 


4.3,2 


[A] is the modal matrix with elements A ^ , where the latin 
subscript refers to the element number and the greek subscript 
to the mode number. The modal matrix [Aj can be multiplied 
by an arbitrary diagonal matrix and the product would still 
satisfy Eqn. 4.3.2. From the symmetry of [k] it follows that 
both [A] T [A] and' [A] T [k] [A] are diagonal matrices . 
Selecting a suitable normalization for [A] one can obtain 


[*f m - m - 

[A] T [k][A] - M 


4.3.3 

4.3.4 
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From 4.3*3: 

W" 1 - M T 


4.3.5 


The inverse modal matrix equals the transposed modal matrix. 
Equations 4.3.2 and 4.3.4 read in subscript notation 


* 2 
V ^ 


= 2 
i 


k ji k iv 


where 




= 


Xl/ 

2 A 




k A 

3* Ji i*' 


4.3.6 

4.3.7 

4.3.8 


is an element ©f the diagonalized stiffness matrix.. 

If all natural frequencies are different, the modal 


columns A 




form a system of linearly, independent vectors 


and the expansion theorem holds, for example 


M = HW ■ or y j -i vv ■ 4.3.9 

The rj v are the modal deflections. 

Inserting. Eqn. 4.3.8 in Eqn. 4.3.1 and premultiplying by M* 
one obtains s 


W T M {*} * [*] T M MW* M T M MW 

[«] T { r } 4.3.10 

Because of Eqn. 4.3.3 and 4.3.4 this reduces to - ' 

M + C a J t E°] MM + P1W ,-W 4 - 3 - u 

where the modal force column {v 5 } is given by 

M= M T {f} ° r v ? v f j 4 * 3,12 

In general the modal equations 4.3.10 are coupled by the cross 
damping terms. Only if [c] is of the form 

[c] » a + b [k] 


4.3.13 
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with arbitrary constants a and b is the triple product in 
Eqn. 4.3.10, because of Eqns, 4.3.3 and 4. 3. 4, a diagonal matrix. 
In this case one can write 


T 


HR 


[‘TMM ■ 

and Eqn. 4.3.10 becomes 


or 


4.3.14 


4.3.15 


%+ 2 *> u r) v = <P U 

If Eqn.' 4.3.12 is not satisfied and if the damping is small, the 

m **. 

off-diagonal terms of the triple product £A ^ [c] [AJ are often 

neglected and the uncoupled equations 4.3.15 are used. 

The modal frequency response function for the i/th mode is 
according to Eqn. 4.3.15 


Hw ( <u) 






so that 


L SJ( ! 


4.3.16 


4.3.17 


Inserting {jj} in Eqn. 4.3.9 and considering Eqn. 4.3.12 one obtains 


I f- i r i ii 

A I Aj , 


;y> = l A | | H v (w) I I a I 

The general admittance matrix for Eqn. 4.3.1 defined by 


4.3.18 


[h( <a )J = £ [l] + io[c] + [k]l'V 

is therefore obtained from the diagonal modal admittance 
matrix [h„ ( *>)] by 

[h(»>] - [»][•*.(. )J [Aj 1 

or H, .( <o ) = £ A. H v ( <u ) A . 


4.3.19 
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For a probabilistic description of the response of 
discrete systems to random excitation, using real normal modes, 
the cross correlations and cross power spectral densities 
must be expressed in terms of normal coordinates* We then have 
to distinguish between the following cross correlation 
functions: (As before, it is assumed that all mean values are 

zero) 


Local 

force correlation 

R Jk<T> * * [ 

>/(t)f k (t + r)J 

4.3.20 

Local 

deflection correlation 

R y k (T) - E [ 

yj*(t)y k (t+r)J 

4.3.21 

Modal 

force correlation 

tC* (T) = E 

*j,(t+r)| 

4.3.22 

Modal 

deflection correlation 

(r) = E 

!7 M (t4T)j 

4.3.23 


In matrix notation these 4 functions are written respectively 
[R f (r)] , [R y (r)j , [r^ ( r)J , Jr* (t)] ( 

Physical forces and deflections are, of course, real quantities. 

However, in the next section a description will be given which 

) ■* . 

admits complex coordinates. For this reason the preceding 
definitions have been formulated in such a way that they also 
apply to complex force and deflection coordinates. 

Inserting eqn. 4.3*9 in eqn. 4.3.21 one obtains with 
eqn. 4.3.23 

py - V A * R* 7 A 

jk % > k* 

T 4.3.24 

- M - [*'] [«”] W 

This is a relation between modal and local cross correlation 
matrix. Eqn. 4.3*24 can be easily inverted because of 
M C A J “ [l]and one obtains 

["'] • [ *’] V] [•] 


4.3.25 
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Similar relations exist' for the force correlation matrices; 

M = 01 [ r ] [ A f H - 3 ‘ 26 

[«*] - [ A *] [ Rf ] [ A ] ' 4 . 3.27 

The cross spectra are obtained from the cross correlations 
by 

s y w > - hr r e ' 1 " rR Jk {r)dT *- 3 - 28 

“CO 

with the inverse 

CO 

R*[ k (0 = 4.3.29 

Therefore, the transformations 4.3*24 and 4*3*25 are also 
valid for the cross spectra: 

M 

H 

Similar relations exist for the force spectra; 

M -M [*] [*] T 

M -[*f FJ [*] 

The modal frequency response function 4.3.16 now relates 
the modal force and modal deflection spectra: 

s’, (-) = H„*(«)S*, (-)H„ («.) 4.3.34 

a relation which is obtained from eqn. 4.2.5 by premultiplication 
with Fa* 1 and post multiplication with [A] and by considering 


[*1 1-1 N 

[*TH M 


4.3.30 

4.3.31 
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eqns . 4.3,19, 4.3.31 and 4.3.33., Rather than manipulating the 
complex admittance matrix [H ( w )] used in eqn. 4.2.5, the 
normal mode analysis allows a very simple computation of res- 
ponse spectra from input spectra according to eqn. 4.3.34. 

If the local force spectra are given and the local deflection 
spectra are desired, it is necessary to first transform the 
local force spectra to modal force spectra with the help of 
eqn. 4.3.33, and then to transform the modal deflection spectra 
to local deflection spectra with the help of eqn. 4.3.30. 

Finally the local cross correlation matrix is calculated 
from eqn. 3.5*8 for t » 0: 

oo 

= 4 . 3.35 

•CO 

The diagonal terms of this matrix are the mean square deflec- , 
tions of the blade elements. It should he once more pointed 
out that the real normal mode analysis is only applicable if 
either the damping matrix has the form given in eqn. 4.3.13 
or if the cross damping terms in the modal equations 4.3.11 

can be neglected. 

4 .4 Complex Normal Mode Analysis 

For blade flap bending the conditions for the validity 
of the real normal mode analysis are not well satisfied, 
since the damping matrix is neither proportional to the mass 
or stiffness matrix nor are the cross damping terms in the 
modal representation negligibly small. A rigorous normal mode 
analysis without neglecting terms is possible if one converts 
the m second order differential equations with respect to time 
into 2m first order equations equivalent to Hamilton’s 
canonical equations. The new system of equations has 2m 
complex eigenvalues and 2n complex eigencolumns, if one 
excludes the case of multiple eigenvalues. It Is further 
assumed, that the system, is stable, though unstable systems 
can also be treated with this method if another coordinate 
transformation is performed, see Ref. (18) . The inconvenience 
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of complex eigencolumns and complex normal coordinates in this 
method is in part compensated by the very simple form of the 
frequency response function. 

With the transformation 



for deflections 4.4.1 



for forces 4.4,2 


equation 4.2.1 assumes the form 

[M] <y| + [k] |y} = jf} 


4.4.3 


where 



4.4.4 


4.4.5 


The eigenvalue problem for eqn. 4.4.3 is defined by 


[M] + K A — 0 4.4.6 


where p] is the modal matrix. Eqn. 4.4.6 is the same as 
eqn. 4.3.1 except for the sign. Prom the symmetry of the 
matrices [M] and [Kj one can derive the orthogonality 
relations 


T 

M 

[M] 

P3 

*7 

II 

T 

m 

[K] 

[5] 

= M 


4.4.8 
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If all eigenvalues \ v are different the complex eigenvectors 
K form a system of linearly independent vectors and the 

j” 

expansion theorem holds: 


7j - E r or-J7f - [*1W 


4.4.9 


Inserting this sum into eqn. 4.4.3 and premultiplying by [5 
one obtains with eqns. 4.4.7 and 4.4.8: 


where 


Mini + MW = M 

{*} = [if |?} or = Efi ir Tj 


4.4.10 

4.4.11 


With 


M I'T-N 


4*4.12 


one can write the set of 2m uncoupled equations 4.4.10 in the 
form 

-1 


’TMM -MW 


or 


4.4.13 


V K ^ 




ft* 


Eqn. 4.4.13 has the modal frequency response function 


H, («) = 


4.4.14 


( lta~kv)ft* 

Formally all relations of section 4.3 are still valid if y 
and Tf are replaced by y and rj and f 1 and f are replaced by 
f an d v respectively. While in section 4.3 the stars on [A] 
were meaningless since (A] was real, the complex character 
of m must now be properly considered. 
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Given the local force spectra 
transforms to modal force spectra by 

[s ? (»)] =[r] [s' {<-)] [ a] , 

then to modal deflection spectra by 

Sjjjj ( w ) » Hjj (<v) ( w ) ( w ) 


[s’ (-)] 


or 


17J 


(-) - 






one first 


4 .A. 15 


4.4.16 


and finally to local deflection spectra by 

T 

■[*.,]- [a*] [s* (.,] [a] 4.4.17 

The transformation of the physical spectra £s y («*»)J and 
[S* («j to the values for the first order system [s ? («)] 
and ¥ < w 0 are to be performed in agreement With equations 
4.4.1 and 4.4.2. 

Since S^ k (w) » E ^Y^(w)Y k (w)J 

one obtains from eqns. 4.4.1 and 4.4.2: 




2 

w 

V(«) ] 

J iw 

} y <“)] 

-iw 

s y (»)] 

T “ “ 

[ sy (“0 



4.4.18 


4.4.19 
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If [s y (w)] is given as the result of the modal analysis, 
it is easy to obtain from eqn. 4.4*18. In Appendix C 

an algorithm is developed for the complex normal mode analysis 
of randomly excited multi-degree of freedom systems, 

. If [S] is the stress-deflection matrix with matrix 


element ( 


<7 k indicating the bending stress at blade 


station k) the physical stress cross correlation matrix is 
given by 


[O ■ M M bl 


4.4.20 
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5. Outline for Lifting Rotor Model Tests In Turbulent Flow 
In the random blade flapping and blade flap bending 
theories presented in Sections 3 and 4, it is assumed that the 
random blade loads at the various blade stations are known, 
so that the output random flapping angles or flap bending 
deflections or moments can be computed. In the following, 
two types of model tests are outlined which have two different 
relations to the theory. In the first and more sophisticated 
type of model -tests, the test results are used to substantiate 
the theory. In the second less sophisticated type of model 
tests the theory Is required to interpret the test results. 

In both types of model tests information is obtained about 
the aerodynamic blade loads when the model Is operated In 
a turbulent flow. 

5.1 Model Tests with Both Load and Load Response Measurements 
In order to substantiate the theory by experiments, 
it is required to record on tape at a number of blade stations 
the aerodynamic blade loads together with flapping angle 
(if flapping hinges are provided) and flap-bending strains. 

The tape recordings must be of sufficient length to allow 
adequate random data processing, say 500 rotor revolutions. 

The pressure pick-ups must be fast responding and the slip 
rings should not Introduce excessive noise. Probably it will 
not be possible to install a sufficient number of pressure 
pick-ups to reliably integrate the blade aerodynamic loads, 
and it will be necessary to calculate the loads from the readings 
of a few strategically located pick-ups, possibly just from 
two pick-ups per blade station installed in the region of 
maximum pressure at the upper and lower airfoil surface. 

This system needs calibration in a flow of known fluctuations 
of the flow direction and of known overall aerodynamic pressures, 
which in Itself will present' problems . In addition to the 
model rotor instrumentation the turbulence characteristics 
of the inflow must be measured upstream of the model rotor, 
for example with a number of crossed hot wire probes. 
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These tests will require considerable preparation, but they 
should be, valuable in obtaining data on both the aerodynamic 
loads and on the load responses of lifting rotors operating 
in turbulent flow, thus making a substantiation of the theory 
of Sections 3 and 4 possible. 

5.2 Model Tests with Load Response Measurements Only 

A much simpler test set-up is obtained, if the aero- 
dynamic load measurements are omitted. In this case, only 
flapping angles and flap-bending strains would be recorded 
at a number of blade stations. The theory would then be used 
in order to establish the random aerodynamic load characteristics 
from the random response measurements. Same as in the previous 
case, the turbulence characteristics of the Inflow must be 
measured upstream of the model rotor. These tests would not 
be a check of the theory which has to be assumed as valid, 
but the tests would be valuable in obtaining data on the random 
aerodynamic blade loads produced by the turbulent inflow. 
Depending on the character of this turbulence it might be 
possible to predict theoretically the rotor blade angle of 
attack fluctuations and thereby the blade load fluctuations 
from a given turbulence spectrum of the inflow. Most likely, 
however, the relation between inflow turbulence and the tur- 
bulence which the rotary wing feels will be a rather complex 
one and may not be easily tractable by a theoretical approach. 

At Ames Research Center there is available a 4-bladed 
12 ft* diameter propeller which could be used as generator 
for a turbulent flow by setting the two blade pairs at largely 
different incidence angles. Also available is a two bladed 
see _ saw rotor of 9 ft. diameter which could be operated in 
the wake of the 4-bladed propeller. 

A Minimum Instrumentation for the test is considered 

as follows: 

1. Flapping angle pick-up, probably best obtained with 
a strain gauged flexure which is loaded in blade 
flapping . 
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2. Two strain gauges near the blade root, one for each 
blade, measuring flap-bending strain. Possibly more 
flap-bending strain gauges along the blade radius. 

3. Blade azimuth indicator. 

4. Two crossed hot wire probes, located upstream of the 
rotor in the rotor plane with a lateral distance of 
about .7 rotor diameters between each other. 

The readings of items 1 to 4 must be tape recorded. In 

addition, blade incidence and rotor shaft incidence must be 


known . 





6 . Conclusions 

g j concepts for a theoretical and experimental study 

of lifting rotor random loads and vibrations have been 
developed. Random lifting rotor loads occur not only during 
operation in atmospheric turbulence but also in a quiet 
atmosphere when there can be a self induced turbulence of the 
flow through the rotor plane as in vertical descent in the 
vortex state/ or as in transition from hovering to forward 
flight. 

6.2 A survey of random process theory has shown that very 

little prior work has been done in the area of non stationary 
random processes, in particular when the system parameters 
are time varying* The only type of non stationary random 
processes which are presently analytically manageable are 
modulated stationary random processes. Under the assumption 
that the random processes associated with lifting rotor 
operation are of this type, it has been shown that up to an 
advance ratio of M = .3 the effect of the time variability 
of the system parameters in blade flapping is quite small. 

If confirmed by experimental hardware or simulator results, 
this would greatly simplify lifting rotor random loads and 
vibrations analyses except for very high advance ratios. 

6 # 3 a rather complete theory of random blade flap-bending 

has been presented under the usual assumptions of stationary and 
ergodic behaviour. While the theoretical concepts are available, 
the numerical analysis is quite extensive. No previous work 
exists which would indicate the validity of possible approxi- 
mations like the neglect of cross damping terms in the modal 
representation, or the neglect of the cross correlations 
between blade elements or blade modes, or the less drastic 
neglect of the quad spectra between blade elements or blade 
modes. A systematic computational and experimental approach 
to these important questions will be required. It is reassuring, 
however, that a complete cross correlation analysis is available 



in principle and may well be computationally tractable. 

Experimental work with lifting rotor models is proposed 
to study the random loads and load responses when the model 
is operating in a turbulent flow. With a rather sophisticated 
rotor model instrumented for the recording of both aerodynamic 
blade loads and blade responses to these loads, the theory 
could be substantiated. However, valuable test results can 
also be expected from a simple rotor model instrumented 
only for the recording of flapping angles and flap -bending 
strains , In such tests the theory would be required to 
interpret the test results with respect to random air loads 
produced by operation of the model rotor in turbulent flow. 
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Appendix A: Time Averaged Autocorrelation Function and 

Associated Power Spectral Density of a 
S tationary Random Process Modulated by a 
Periodic Function 

if y(t) is the product of a stationary random function 
x(t) and a deterministic periodic function A(t) (modulated 
stationary random process), then the random process obtained 
by averaging the autocorrelation function 

c* 

T ^Tl R y( t > T ) dt = R y < T > *-1 

~t/ 2 

is a stationary random process and has the measurable single 
frequency power spectral density, see Section 9,9.3 of Kef. O). 

a> 

s y (w) = j R y (r)e“ iwT dr A -2 

—a) 

As A(t) is periodic, it can be expressed in Fourier series 

A(t) = £ 

k =-® 

Setting 

y(t) = x(t) £ c k^ k ‘ b ' 

k 

one obtains 

Ry ( ¥2 )=E[x(t,)x(t 2 ) { £ 

k eH» |= — CO J 

k b»0D I =— CD 

tl+tp 

Substitute = r and — ^ — = t, to get 

ynt) = ^(r) £ g Ck0| .M(k+lrt+(l-kni» 

k I «-® 
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While applying equation A-l, it is easy to verify that 


lim I 
T-~ ® T 


_ 


9 

i(il 0 (k-f f ) 




except for k = I 
& k * ti 



As I and k can take both positive and negative values I -k 
can take values 2k and ~2k. 

Therefore, 


B y (r) - 


■Ur) S =/ 

A —co 


°o R x< T > 


k * 0 

Applying equation A-2, one gets a measurable single frequency 
power spectral density 


s y («) » S x (w)c 0 2 + £ c k 2 | S x (w + S x («-k60l A-3 

k *- o 

Another way of obtaining the same result is to multiply the 
double frequency spectrum cor reponding to the sample function 
x(t)A(t) by i(w 2 -w, ) . This is shown in the following deriva- 
tion. 

With 

y(t) - x(t) 

k =-oo 

set 

Y*(w,) = j x( t , ) 

*" CD k =-co 

to obtain 

Y*(w 1 ) =c 0 X 4 '(w 1 ) + C k * [ x*(w 1 -kw 0 ) + X*l[w^+kfe> e )J 

— oo 

k * o 
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and 


Y( 0> 2 ) = c Q X( w 2 )+ E ° k [ x ( w 2- k « 0 ) +x ( w 2 +k "p] 


— 00 

k t 0 


The double frequency spectrum is defined by 

Sy( , 0> 2 ) “ E [y ( W^)Y( W 2 )J 

We now multiply this double frequency spectrum with $( Uj ) 

2 Cilg 

«( W 2 - wpSyC w 2 ) ■= *( W 2 - S x ( — ? — ) 


cc w p -2ko»+w, 

51 o k «( w 2 " W 1^ S x^" 2 * 


— 00 

k * 0 
00 


A -4 


+ E c k 2 *( W fT W 1 )S x ( 
■ —00 


«i# 2 +2k fcl 0 + 


k * 0 


and obtain the same equation as A-3 
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Appendix B: Computation of Power Spectr al Density f or _a 

Flapping Randomly Excited Blade 
Equation B-1 represents a general type of linear, 
second order and non-homogeneous differential equation with 
constant and periodic coefficients, 0(t) and ®(t) 
representing respectively the non stationary random response 
and the stationary random input. 

0 + { c-j+a-j sin( nw 0 t +6^)} 0 + { c 2 +a 2 Sin ^ n<4 ^ 

*= 1 c 3+ a 3 sin ( n H f + e 3 ) } Q! ( t ) 

(Selecting a time unit for which the rotor angular velocity^ 

nw,= 1 and substituting ~ 0, © 2 = = 0, = c 3 g> 

a „ a ■ c 2 « 1, and a 3 - ^ in equation B-1, an approximate 

blade flapping differential equation valid up to about **=0.3 
advance ratio, is obtained. For details refer to Section 3.) 

With the operator f representing the Fourier integral trans- 
formation, let 

f( $) «. - o> 2 b( w ) 

f|jc|+a.|Sln(n%t +©-j )}^j .« ( w ) 

^[{c 2 +a 2 sin( no^t +© 2 ^] = B 2 ( « ) and 

^[{ c 3 +a 3 sin ( n ^ m B 3^ " ) 

Apply the operator f to equation B-1 at frequencies <*] and 
to yield 

-t ^ 2 )+B 2 *(w 1 ) ® B 3 (w-j) 

and -w 2 2 B(w 2 )+B 1 (w 2 )+B 2 (w 2 ) » B 3 (<i* 2 ) . 
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Now, consider 

e[b*( • 1 )B( « 2 )j 

along the line ^ — w 2 “ u 9 as Appendix A 

and set 

$( w 2 " wj)e £b*( w^)B( = % (**')• 

*1 

After some algebra one obtains 

|w^+Cj^ w^- 2 Cg w ^ +c 2 ^ } S# ( w ) + { 4 ” ( w “ nu *) " 2 ( w -n«o ) sin( 9<j “&>) 

a ^ . 

+ Tjj 1 " j ( ^—nw^ ) 

2 a 2 
+ i !^_( u+nn ,) 2 + «+nw 0 )sin(e 1 -e 2 )+ ( « +n«b ) 

a 2 ' 

= o^S a ( w ) 4 - 4 " | S a ( w +nw 0 )+S a ( w -nw© ) | B-2 

With a> = ko > 0 where k is a discrete variable with values k *= 1,2,... 
and 6 >o a fixed step size, eqn. B-2 represents an infinite set 
of difference equations for (k<u 0 ) which will be abbreviated 
by S k < w ).. Since S k _ n ( » ) = S. k+n ( ) and since from physical - 
considerations S k+n ( ) — ** 0 for large k, the infinite set of 
equations can be approximated by a finite set. 

Express equation B-2 iri a symbolic form 

A(k,k)S k ( <»> )+A(k-n,k)S k n ( a) )+A(k+n,k)S k+n ( w ) « R.H.S. B-3 

Equation B-3 is inverted for S k ( cj ) values by a digital computer 
program ( 19 ) which essentially involves a matrix inversion 
subroutine. Figures 2o-d show numerical results for two stochastic 
Inputs ( g>) . Figures 2a and 2b with a^ = a 2 .» .2 refer to 
an advance ratio n = 0*3. and a blade inertia number Y *= 4. 

For this set of coefficients the values A(k-n,k) and A(k+n,k) 
are two order of magnitude smaller than A(k,k) which explains 
the fact that the results of the present algorithm do not 
differ appreciably from those for the constant parameter system 
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obtained by omitting the periodic terms in eqn. B-1 . The 
results obtained for the same set of stochastic inputs but , 
now with a-j — a^ — 0*6 show a substantial difference* 

Figures 2c and 2d, between the present method and the equivalent 
constant parameter system, as expected. It is to be noted 
that the second set of coefficients have no physical relevance 
as the approximate blade flapping equation is valid up to 
about = 0.3 advance ratio and the blade inertia number 
is generally less than 10 for helicopters. 




FIGURE 2,: SPECTRAL DENSITY OF FLAPPING OSCILLATIONS OF RIGID SLADES IN FORWARD 
HELICOPTER FLIGHT, ASSUMING THE RESPONSE TO BE A STATIONARY RANDOM 
PROCESS MODULATED BY PERODIC FUNCTIONS. (Equations given above 
represent the corresponding blade flapping equation and the assumed 
spectral density of the mean blade angle of attack 3 respectively.) 
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Appendix C : Discrete Element Algorithm for Complex Normal 

Mode Analysis of Multi Degree of Freedom Stable 
Systems Under Random Excitation 
Before actually coming to the numerical analysis,, it 
is mandatory to replace the continuum as an assemblage of 'm» 
discrete elements to obtain ’m* equations of dynamic equilibrium 
in a standard form, equation 4, .2.1, and then, as discussed in 
Section 4.4, convert these m equations into 2m first order 
differential equations having 2m complex eigenvalues and 2m 
complex modal columns. The discussion here is restricted 
to blade flap-bending random analysis with small advance ratio 
though the present algorithm is applicable to any discrete 
systems under random excitation. 

The analysis of numerical problems comprises three 

main phases: 

1) the mass, damping and stiffness matrices 

2) the eigenvalues and eigencolumns 

3) the spectral density matrices and RMS deflection 

and stress values at element stations. The numerical 
operations are carried out in two passes - the first 
two phases In the first run and the last one in the 
second run . 

The computer program, written in FORTRAN IV language 
is based on some of the subroutines appearing elsewhere in 
the literature (19, 20, 21, 22) . It accepts as inputs the 
following information regarding the pattern of discretization 
and other constants associated with the physical system. 

i) number of elements 

ii) Hubradius 

iii) boundary conditions at the hub joint 

iv) lift slope (5.7 per radian) 

v) blade chord at midpoint of the element or element 
station 

vi) air density 

vii) rotor speed in rad/sec 

viii) gravitational constant 

ix) X-ordinate of element stations 
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x) flexured rigidity El at element stations 

xi) vertical load at element stations 

In the present case, the mass and damping matrices 

are diagonal matrices, hence they can be directly read into 

the program. The damping terms at each element station can 

also be calculated from the formula ac p ( see equation 4.1.1). 

Depending upon the initial conditions at the hubjolnt solve for 

Mi, * — M . and 8. ,0 O , from equations 

1 * 2 m- 1 i'2 m 

4.1.6 and 4.1.8 by calling a matrix subroutine which computes 
the inverse and rank of a matrix and solves a set of simultaneous 
equations according to Gauss-Jordan reduction principle. 

Generate the matrix of influence coefficients, W^, from the 
recurrence relation 4.1.9. The stiffness matrix^k^jJ is 

obtained by inverting the matrix M- ..... 

With Qn] , [cj , and DO ‘matrices being generated in 
the first operation, call for matrix subroutines to compute 


L-I'M 



and finally compute the eigenvalue matrix (see equation 3.6.8) 






[•rw 

[*: 


[» 



The 2m eigenvalues and 2m eigenvectors are calculated in the 
following order: 

i) Generating the characteristic polynomlnal. The subroutine 

is based on the Leverrier-Faddev method (20-355) 
according to which the coefficients P| ,pg . • ,P n in the 
characteristic polynomlnal 


v n „ 

X "P i 



n-2 

X 



0 
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are given by 

P n = (l/n)( trace ) 

where 

[«„]*[«] [K-a-'.-. W] 

il) Calculation of roots of the characteristic polynominal: 

This is effected by a subroutine of reference 20 
which finds real and imaginary parts of complex roots 
and gives the reduced polynominal vjhich in turn is 
examined for real roots if any. 

iii) Solution of simultaneous linear equations for eigenvectors f 
The complex modal columns or eigenvectors are obtained 
by a simultaneous equation subroutine ( 19 > 22). 

Presuming that the physical cross spectra IS ■ ( «*> )J. 

of the exciting system are known, the RMS deflection and stress 

values are computed in the following sequence: 

1) Obtain the local force cross spectra from equation 4 . 4.19 
and then transform it to the modal force cross spectra 
according to the relation 4 . 4 . 15 . 

2) Prom equation 4,4.16 compute the modal deflection cross 
spectra which when substituted in equation 4.4.17 
gives the local deflection cross spectra. 

3) Obtain the physical deflection cross spectra from 
equation 4.4.18. 

4) Using any quadrature subroutine calculate the physical 
deflection and stress space cross correlation 

[R y (0)j - J [s y («)] dw 
[ r ° ( 0 )] « [s] [R y (0)] [s] 

the diagonal terms of which give the corresponding RMS 
values at element stations (refer to equations 4.4.20). 



